!----------------------------------------------------------------------------
! Program: 2D_Potential_FMM - A fast multipole boundary element
! Method (BEM) code for analyzing large-scale, general 2D
! potential problems (governed by Laplace equation)
! using constant elements.
!
! Developers: Dr. Naoshi Nishimura at Kyoto University, Japan;
! Dr. Yijun Liu at the University of Cincinnati, Cincinnati, OH.
!
! Version: V.1.20.
! Released: October 1, 2008.
!
! Copyright(c)2004-2008 By Kyoto University and University of Cincinnati.
! This code is intended for educational use only. No part of
! the code can be used for any commercial applications/
! distributions without prior written permissions of the
! original developers.
!
!----------------------------------------------------------------------------

program fmm_2d_potential
   implicit real*8(a-h, o-z)
   integer, allocatable:: ia(:)
   complex*16, allocatable :: am(:)
   character*80 Prob_Title
   call CPU_Time(time0)

   open (4, file='input.fmm',    status='old')
   open (5, file='input.dat',    status='old')
   open (3, file='output.dat',   status='unknown')
   open (7, file='phi_boundary.plt', status='unknown')
   open (8, file='xy.plt',       status='unknown')
   open (9, file='phi_domain.plt',     status='unknown')
   open (10, file='debug.dat',   status='unknown')

! Input the parameters
   read(4,*)   maxl, levmx, nexp, ntylr, tolerance
   read(4,*)   maxia, ncellmx, nleafmx, mxl, nwksz
   read(5,'(a80)')   Prob_Title
   read(5,*)   n, nfield
   write(3,'(a80)') Prob_Title
   write(*,'(a80)') Prob_Title

! Estimate the maximum numbers of the cells and leaves,
! and size of the preconditioning matrix, etc.
   if(ncellmx.le.0) ncellmx = max(4*n/maxl,100)
   if(nleafmx.le.0) nleafmx = max(ncellmx/2,100)
   if(nwksz.le.0) nwksz = maxl*maxl*nleafmx
   ligw = 20
   lrgw = 1+n*(mxl+6)+mxl*(mxl+3)
   iwksz = n+3*nleafmx+1
   allocate (ia(maxia))

! Load the addresses (pointers) associated with the locations of the
! variables to be stored in the large array "am"
   call lpointer(lp, ln, maxia, ia, n, nexp, ntylr, ncellmx,&
      levmx, ligw, lrgw, nwksz, iwksz, nfield,&
      l_n, l_x, l_y, l_node, l_dnorm,&
      l_bc, l_a, l_b, l_xmax,&
      l_xmin, l_ymax, l_ymin, l_ielem, l_itree,&
      l_level, l_loct, l_numt, l_ifath, l_lowlev,&
      l_maxl, l_levmx, l_nexp, l_ntylr, l_tolerance,&
      l_ncellmx, l_nleafmx, l_mxl, l_u, l_ax,&
      l_sb, l_sx, l_ligw, l_lrgw, l_igwk,&
      l_rgwk, l_nwksz, l_iwksz, l_rwork, l_iwork,&
      l_xfield, l_nfield, l_f)

! Estimate the memory usage

   maxa = lp
   write(3,100) maxa*16/1.D6
   write(*,100) maxa*16/1.D6
100 format(/'Memory size of the large block am =', f12.1,' Mb'/)

! Allocate the large block 'am'
   allocate (am(maxa))

! Assign the parameters to the array am()
   call assigni(n, am(l_n))
   call assigni(maxl, am(l_maxl))
   call assigni(levmx, am(l_levmx))
   call assigni(nexp, am(l_nexp))
   call assigni(ntylr, am(l_ntylr))
   call assignd(tolerance, am(l_tolerance))
   call assigni(ncellmx, am(l_ncellmx))
   call assigni(nleafmx, am(l_nleafmx))
   call assigni(mxl, am(l_mxl))
   call assigni(ligw, am(l_ligw))
   call assigni(lrgw, am(l_lrgw))
   call assigni(nwksz, am(l_nwksz))
   call assigni(iwksz, am(l_iwksz))
   call assigni(nfield, am(l_nfield))

! Call the FMM BEM main program
   call fmmmain(maxa, maxia, am, ia,&
      am(l_n), am(l_x), am(l_y), am(l_node),&
      am(l_dnorm), am(l_bc), am(l_a), am(l_b),&
      am(l_xmax), am(l_xmin), am(l_ymax), am(l_ymin),&
      am(l_ielem), am(l_itree), am(l_level), am(l_loct),&
      am(l_numt), am(l_ifath), am(l_lowlev), am(l_maxl),&
      am(l_levmx), am(l_nexp), am(l_ntylr), am(l_tolerance),&
      am(l_ncellmx), am(l_nleafmx), am(l_mxl), am(l_u),&
      am(l_ax), am(l_nfield), am(l_xfield), am(l_f),&
      am(l_sb), am(l_sx), am(l_igwk), am(l_rgwk),&
      am(l_ligw), am(l_lrgw), am(l_nwksz), am(l_iwksz),&
      am(l_rwork), am(l_iwork))

! Estimate the total CPU time
   call CPU_Time(time)
   write(3,*)
   write(*,*)
   write(3,*) ' Total CPU time used =', time-time0, '(sec)'
   write(*,*) ' Total CPU time used =', time-time0, '(sec)'
   stop
end program fmm_2d_potential

!----------------------------------------------------------------------------
! Definition of Variables:
!
! maxa = maximum size of the array am
! maxia = maximum number of variables allowed
! am = a large array storing the variables for the SLATEC GMRES solver
! ia = an array storing the locations of the variables in the array am
!
! n = number of elements (= number of nodes)
! x = coordinates of the nodes
! y = coordinates of the end points of the elements
! node = element connectivity
! dnorm = normal at each node
! bc = BC type and value
!
! a = multipole expansion moments
! b = local expansion coefficients
! xmax,xmin = maximum and minimum x coordinate
! ymax,ymin = maximum and minimum y coordinate
! ielem = ielem(i) gives the original element number for i-th element in
! the quad-tree structure c itree = itree(c) gives the cell location of c-th cell within each
! tree level
! loct = elements included in the c-th cell are listed starting at
! the loct(c)-th place in the array ielem
! numt = numt(c) gives the number of elements included in the c-th cell
! ifath = ifath(c) gives the number of the parent cell of the c-th cell
! level = level l cells start at the level(l)-th cell in the tree
! lowlev = number of the tree levels
!
! maxl = maximum number of elements allowed in a leaf
! levmx = maximum number of levels allowed in the tree structure
! nexp = number of terms in multipole expansion
! ntylr = number of terms in local expansion
! tolerance = GMRES solution convergence tolerance
! ncellmx = maximum number of cells allowed in the tree
! nleafmx = maximum number of leaves allowed in the tree
! mxl = maximum dimension of Krylov subspace (used in GMRES)
!
! u = first stores b vector; then solution vector of system Ax = b
! ax = resulting vector of multiplication Ax
! nfield = number of the field points inside the domain
! xfield = coordinates of the field points inside the domain
! f = values of the potential at the field points inside the domain
!
! The following variables and arrays are used in the SLATEC GMRES solver:
! sb,sx,igwk,rgwk,ligw,lrgw,nwksz,iwksz,rwork,iwork
!
!----------------------------------------------------------------------------

subroutine lpointer(lp, ln, maxia, ia, n, nexp, ntylr, ncellmx,&
   levmx, ligw, lrgw, nwksz, iwksz, nfield,&
   l_n, l_x, l_y, l_node, l_dnorm,&
   l_bc, l_a, l_b, l_xmax,&
   l_xmin, l_ymax, l_ymin, l_ielem, l_itree,&
   l_level, l_loct, l_numt, l_ifath, l_lowlev,&
   l_maxl, l_levmx, l_nexp, l_ntylr, l_tolerance,&
   l_ncellmx, l_nleafmx, l_mxl, l_u, l_ax,&
   l_sb, l_sx, l_ligw, l_lrgw, l_igwk,&
   l_rgwk, l_nwksz, l_iwksz, l_rwork, l_iwork,&
   l_xfield, l_nfield, l_f)

   dimension ia(maxia)
   lp= 1

   l_n = l_address( 1,maxia,ia,lp,4,1)
   l_x = l_address( 2,maxia,ia,lp,8,n*2)
   l_y = l_address( 3,maxia,ia,lp,8,n*2)
   l_node = l_address( 4,maxia,ia,lp,4,n*2)
   l_dnorm = l_address( 5,maxia,ia,lp,8,n*2)
   l_bc = l_address( 6,maxia,ia,lp,8,n*2)
   l_a = l_address( 7,maxia,ia,lp,16,(nexp+1)*ncellmx)
   l_b = l_address( 8,maxia,ia,lp,16,(ntylr+1)*ncellmx)
   l_xmax = l_address( 9,maxia,ia,lp,8,1)
   l_xmin = l_address(10,maxia,ia,lp,8,1)
   l_ymax = l_address(11,maxia,ia,lp,8,1)
   l_ymin = l_address(12,maxia,ia,lp,8,1)
   l_ielem = l_address(13,maxia,ia,lp,4,n)
   l_itree = l_address(14,maxia,ia,lp,4,ncellmx)
   l_level = l_address(15,maxia,ia,lp,4,levmx+1)
   l_loct = l_address(16,maxia,ia,lp,4,ncellmx)
   l_numt = l_address(17,maxia,ia,lp,4,ncellmx)
   l_ifath = l_address(18,maxia,ia,lp,4,ncellmx)
   l_lowlev = l_address(19,maxia,ia,lp,4,1)
   l_maxl = l_address(20,maxia,ia,lp,4,1)
   l_levmx = l_address(21,maxia,ia,lp,4,1)
   l_nexp = l_address(22,maxia,ia,lp,4,1)
   l_ntylr = l_address(23,maxia,ia,lp,4,1)
   l_tolerance = l_address(24,maxia,ia,lp,8,1)
   l_ncellmx = l_address(25,maxia,ia,lp,4,1)
   l_nleafmx = l_address(26,maxia,ia,lp,4,1)
   l_mxl = l_address(27,maxia,ia,lp,4,1)
   l_u = l_address(28,maxia,ia,lp,8,n)
   l_ax = l_address(29,maxia,ia,lp,8,n)
   l_sb = l_address(30,maxia,ia,lp,8,n)
   l_sx = l_address(31,maxia,ia,lp,8,n)
   l_ligw = l_address(32,maxia,ia,lp,4,1)
   l_lrgw = l_address(33,maxia,ia,lp,4,1)
   l_igwk = l_address(34,maxia,ia,lp,4,ligw)
   l_rgwk = l_address(35,maxia,ia,lp,8,lrgw)

   l_nwksz = l_address(36,maxia,ia,lp,4,1)
   l_iwksz = l_address(37,maxia,ia,lp,4,1)
   l_rwork = l_address(38,maxia,ia,lp,8,nwksz)
   l_iwork = l_address(39,maxia,ia,lp,4,iwksz)
   l_xfield = l_address(40,maxia,ia,lp,8,nfield*2)
   l_nfield = l_address(41,maxia,ia,lp,4,1)
   l_f = l_address(42,maxia,ia,lp,8,nfield)
   return
end

!----------------------------------------------------------------------------

integer function l_address(ln,maxia,ia,lp,ibyte,length)
   dimension ia(maxia)
   l_address = lp
   ia(ln) = lp

   iu = 16
   inc = (ibyte*length-1)/iu+1
   lp = lp+inc
   if(ln .gt. maxia) then
      write(*,*)'!Specified # of variables maxia',maxia,'is too small'
      stop
   endif
   return
end

!----------------------------------------------------------------------------
subroutine assigni(i,ii)
   integer i,ii
   ii = i
   return
end

subroutine assignd(d,dd)
   real*8 d,dd
   dd = d
   return
end

!----------------------------------------------------------------------------
subroutine fmmmain(maxa, maxia, am, ia, n,x,y,node,dnorm,bc,&
   a,b, xmax,xmin,ymax,ymin,ielem,itree,level,loct,numt,&
   ifath,lowlev,maxl,levmx,nexp,ntylr,tolerance,ncellmx,&
   nleafmx,mxl,u,ax,nfield,xfield,f,sb,sx,igwk,rgwk,&
   ligw,lrgw,nwksz,iwksz,rwork,iwork)

   implicit real*8(a-h,o-z)
   complex*16 am(maxa), a,b
   dimension ia(maxia),ja(1),a(0:nexp,ncellmx),b(0:ntylr,ncellmx),&
      x(2,n),y(2,n),node(2,n),dnorm(2,n),bc(2,n),&
      ielem(n),itree(ncellmx),level(0:levmx),loct(ncellmx),&
      numt(ncellmx),ifath(ncellmx), u(n),ax(n),sb(n),sx(n),&
      igwk(ligw),rgwk(lrgw),rwork(nwksz),iwork(iwksz),&
      xfield(2,nfield),f(nfield)

   external matvec, msolve

! Input parameters and prepare the BEM model
   call prep_model(n,x,y,node,bc,dnorm,xfield,nfield,maxl,levmx,&
      nexp,ntylr,tolerance,xmin,xmax,ymin,ymax)

! Generate the quad-tree structure for the elements
   call tree(n,x,xmax,xmin,ymax,ymin,ielem,itree,level,loct,numt,&
      ifath,lowlev,maxl,levmx,ncellmx,nleafmx,nwksz,iwork)

! Compute the right-hand-side vector b with the FMM
   call fmmbvector(n,x,y,node,dnorm,bc,u,ax,a,b,xmax,xmin,ymax,ymin,&
      ielem,itree,level,loct,numt,ifath,&
      nexp,ntylr,ncellmx,lowlev,maxl,rwork,iwork)

!    write (10,*) "moment arr a"
!    write(10,101) ((a(i,j),j=1,nexp+1),i=1,iwork(1))
! 101 format(15(ES13.6,2x))

!    write (10,*) "L coes arr b"
!    write(10,102) ((b(i,j),j=1,ntylr+1),i=1,iwork(1))
! 102 format(15(ES13.6,2x))

!    write (10,*) "moment arr a"
!    write (10,*) a
!    write (10,*) "ls b"
!    write (10,104) u
! 104 format(20(ES13.6,2x))
! Solve the BEM system of equations Ax=b with the fast multipole BEM

! Prepare parameters for calling the iterative solver GMRES
! (SLATEC GMRES solver is used, which is available at www.netlib.org.
! See the documentation for the SLATEC GMRES solver for more information
! about the following related parameters)

   nelt = 1
   isym = 0
   itol = 0
   tol = tolerance
   iunit = 3
   igwk(1) = mxl
   igwk(2) = mxl
   igwk(3) = 0
   igwk(4) = 1
   igwk(5) = 10
   do i=1,n
      ax(i) = 0.d0
   enddo

   write(*,*) 'Call Equation Solver GMRES ...'
   call dgmres(n,u,ax, nelt,ia,ja,am,isym, matvec,msolve,itol,tol,&
      itmax,iter,er,ierr,iunit,sb,sx,rgwk,lrgw,igwk,ligw,&
      rwork,iwork)
   write(3,*) ' Error indicator from GMRES:', ierr
   write(*,*) ' Error indicator from GMRES:', ierr

! Output the boundary solution
   do i=1,n
      u(ielem(i)) = ax(i)
   enddo

   write(3,*) ' Fast Multipole BEM Solution:'
   do i=1,n
      write(3,*) i, u(i)
      write(7, '(I9, 4X, ES13.6)') i,u(i)
   enddo

! Evaluate the field inside the domain and output the results
   call domain_field(nfield,xfield,f,n,x,y,bc,node,dnorm,u)
   return
end
!----------------------------------------------------------------------------

subroutine prep_model(n,x,y,node,bc,dnorm,xfield,nfield,maxl,&
   levmx,nexp,ntylr,tolerance,&
   xmin,xmax,ymin,ymax)

   implicit real*8(a-h,o-z)
   dimension x(2,*),y(2,*),node(2,*),bc(2,*),dnorm(2,*),xfield(2,*)
   write(*,2) n, maxl, levmx, nexp, tolerance
   write(3,2) n, maxl, levmx, nexp, tolerance
2  format(' Total number of elements =', I12&
      /' Max. number of elements in a leaf =', I12&
      /' Max. number of tree levels =', I12&
      /' Number of terms used in expansions =', I12&
      /' Tolerance for convergence =', D12.3)
   write(*,*)
   write(3,*)

! Input the mesh data
   read(5,*)
   do i=1,n
      read(5,*) itemp, y(1,i), y(2,i)
   enddo
   read(5,*)
   do i=1,n
      read(5,*) itemp, node(1,i), node(2,i), bc(1,i), bc(2,i)
   enddo

! Input the field points inside the domain
   if (nfield .gt. 0) then
      read(5,*)
      do i=1,nfield
         read(5,*) itemp, xfield(1,i), xfield(2,i)
      enddo
   endif

! Compute mid-nodes and normals of the elements
   do i=1,n
      x(1,i) = (y(1,node(1,i)) + y(1,node(2,i)))*0.5
      x(2,i) = (y(2,node(1,i)) + y(2,node(2,i)))*0.5
      h1 = y(2,node(2,i)) - y(2,node(1,i))
      h2 = -y(1,node(2,i)) + y(1,node(1,i))
      el = sqrt(h1**2 + h2**2)
      dnorm(1,i) = h1/el
      dnorm(2,i) = h2/el
   enddo

! Determine the square bounding the problem domain (Largest cell used in FMM)
   xmin=x(1,1)
   xmax=x(1,1)
   ymin=x(2,1)
   ymax=x(2,1)

   do 10 i=2,n
      if(x(1,i).le.xmin) then
         xmin=x(1,i)
      elseif(x(1,i).ge.xmax) then
         xmax=x(1,i)
      endif
      if(x(2,i).le.ymin) then
         ymin=x(2,i)
      elseif(x(2,i).ge.ymax) then
         ymax=x(2,i)
      endif
10 continue

   scale = 1.05d0 ! Make the square slightly larger
   xyd = max(xmax-xmin,ymax-ymin)/2.d0
   xyd = xyd*scale
   cx = (xmin+xmax)/2.d0
   cy = (ymin+ymax)/2.d0
   xmin = cx-xyd
   xmax = cx+xyd
   ymin = cy-xyd
   ymax = cy+xyd

! Output nodal coordinates for plotting
   do i = 1,n
      write(8,*) x(1,i), x(2,i)
   enddo
   return
end

!----------------------------------------------------------------------------
subroutine tree(n,x,xmax,xmin,ymax,ymin,ielem,itree,level,loct,&
   numt,ifath,lowlev,maxl,levmx,ncellmx,nleafmx,&
   nwksz,iwork)

   implicit real*8(a-h,o-z)
   complex*16 a,b
   dimension x(2,*),ielem(*),itree(*),level(0:*),loct(*),numt(*),&
      ifath(*), iwork(*), nwk(4)

   do i=1,n
      ielem(i) = i ! Store the original element numbers in ielem
   enddo

! For the level 0 cell (largest cell)
   itree(1) = 0
   level(0) = 1
   level(1) = 2
   loct(1) = 1
   ifath(1) = 1
   numt(1) = n
   ndivx = 1
   lowlev = 1
   nleaf = 0
   nswa = 0

! For cells on level 1 to the lowest level (leaves)
   do 10 lev=1,levmx
      levp = lev-1
      levn = lev+1
      level(levn) = level(lev)
      if(level(lev).eq.level(levp)) goto 900
      ndivxp = ndivx
      ndivx = 2*ndivxp

      dxp = (xmax-xmin)/ndivxp ! Parent cell size
      dyp = (ymax-ymin)/ndivxp
      do 11 inp=level(levp),level(lev)-1
         itrp = itree(inp)
         if(numt(inp).gt.maxl .or.&
            (lev.le.2 .and. numt(inp).ne.0) ) then
            itrpx = mod(itrp,ndivxp)
            itrpy = itrp/ndivxp
            xsep = xmin+(itrpx + 0.5d0)*dxp
            ysep = ymin+(itrpy + 0.5d0)*dyp

            call bisec(x,ielem(loct(inp)),numt(inp),ysep,nsepy, 2)
            call bisec(x,ielem(loct(inp)),nsepy-1,xsep,nsepx1,1)
            call bisec(x,ielem(loct(inp)+nsepy-1),&
               numt(inp)-nsepy+1,xsep,nsepx2,1)

            nwk(1) = nsepx1-1
            nwk(2) = nsepy-nsepx1
            nwk(3) = nsepx2-1
            nwk(4) = numt(inp)-nsepy-nsepx2+2
            locc = loct(inp)

            do 12 icldy=0,1
               do 12 icldx=0,1
                  icld = icldy*2+icldx+1

                  if(nwk(icld).gt.0) then
                     nrel = level(levn)
                     if(nrel.gt.ncellmx) then
                        write(*,*) " ncellmx error"
                        stop
                     endif
                     itree(nrel) = ((itrpy*2+icldy)*ndivxp + itrpx)*2+icldx

                     loct(nrel) = locc
                     numt(nrel) = nwk(icld)
                     ifath(nrel) = inp
                     lowlev=lev

! Leaves:

                     if((lev.ne.1) .and.&
                        (numt(nrel).le.maxl .or. lev.eq.levmx)) then
                        nleaf = nleaf+1
                        if(nleaf.gt.nleafmx) then
                           write(*,*) " nleafmx error"
                           stop
                        endif
                        nleaf3 = nleaf*3 - 1
                        iwork(nleaf3) = nrel ! Store cell number (icell)
                        iwork(nleaf3+1) = nswa + 1 ! Location of pre-cond'er
                        iwork(nleaf3+2) = 1 ! Initial value of switch isw
                        nswa = nswa + numt(nrel)**2
                        if(nswa.gt.nwksz) then
                           write(*,*) " nwksz error"
                           stop
                        endif
                     endif

                     level(levn) = nrel + 1
                     locc = locc + nwk(icld)
                  endif
12          continue
         endif
11    continue
10 continue

900 iwork(1) = nleaf ! Store number of leaves in iwork(1)
   write(3,*) ' Number of tree levels =', lowlev
   write(*,*) ' Number of tree levels =', lowlev
   write(3,*) ' Number of leaves =', nleaf
   write(*,*) ' Number of leaves =', nleaf
   write(3,*) ' Number of cells =', nrel
   write(*,*) ' Number of cells =', nrel
   write(3,*)
   write(*,*)
   return
end

!----------------------------------------------------------------------------
subroutine bisec(x,ielem,n,xsep,nsep,ic)
   implicit real*8(a-h,o-z)
   dimension x(2,*),ielem(*)
   nsep = 1
   if(n.le.0) return

   do ifr=1,n
      if(x(ic,ielem(ifr)).le.xsep) then
         if(ifr.ne.nsep) then
            istore = ielem(nsep)
            ielem(nsep) = ielem(ifr)
            ielem(ifr) = istore
         endif
         nsep = nsep + 1
      endif
   enddo

   return
end

!----------------------------------------------------------------------------
subroutine fmmbvector(n,x,y,node,dnorm,bc,u,ax,a,b,xmax,xmin,&
   ymax,ymin,ielem,itree,level,loct,numt,ifath,&
   nexp,ntylr,ncellmx,lowlev,maxl,rwork,iwork)
   implicit real*8(a-h,o-z)
   complex*16 a,b
   dimension a(0:nexp,ncellmx),b(0:ntylr,ncellmx),&
      x(2,*),y(2,*),node(2,*),dnorm(2,*),bc(2,*),u(*),ax(*),&
      ielem(*),itree(*),level(0:*),loct(*),numt(*),ifath(*),&
      rwork(*),iwork(*)

! Switch the BC type
   do i=1,n
      if(bc(1,i) .eq. 1.) then
         bc(1,i) = 2.d0
      else
         bc(1,i) = 1.d0
      endif
   enddo

   do i=1,n
      u(i)= bc(2,ielem(i))
      ax(i) = 0.d0
   enddo

! Apply the FMM to compute the right-hand side vector b

   call upward(u,n,y,node,dnorm,bc,a,xmax,xmin,ymax,ymin,ielem,&
      itree,level,loct,numt,ifath,nexp,ncellmx,lowlev,maxl)

   call dwnwrd(u,ax,n,x,y,node,dnorm,bc,a,b,xmax,xmin,ymax,ymin,&
      ielem,itree,level,loct,numt,ifath,nexp,ntylr,ncellmx,&
      lowlev,maxl,rwork,iwork)

! Store b vector in u and switch the BC type back
   do i=1,n
      u(i) = - ax(i)
      if(bc(1,i) .eq. 1.) then
         bc(1,i) = 2.d0
      else
         bc(1,i) = 1.d0
      endif
   enddo
   return
end

!----------------------------------------------------------------------------
subroutine matvec(n,u,ax,nelt,ia,ja,am,isym)
   implicit real*8(a-h,o-z)
   complex*16 am, a, b
   dimension u(*),ax(*),ia(*),ja(*),am(*)

! Retrieve the pointers
   l_n = ia(1)
   l_x = ia(2)
   l_y = ia(3)
   l_node = ia(4)
   l_dnorm = ia(5)
   l_bc = ia(6)
   l_a = ia(7)
   l_b = ia(8)
   l_xmax = ia(9)
   l_xmin = ia(10)
   l_ymax = ia(11)
   l_ymin = ia(12)
   l_ielem = ia(13)
   l_itree = ia(14)
   l_level = ia(15)
   l_loct = ia(16)
   l_numt = ia(17)
   l_ifath = ia(18)
   l_lowlev = ia(19)
   l_maxl = ia(20)
   l_levmx = ia(21)
   l_nexp = ia(22)
   l_ntylr = ia(23)
   l_tolerance = ia(24)
   l_ncellmx = ia(25)
   l_nleafmx = ia(26)
   l_mxl = ia(27)
   l_u = ia(28)
   l_ax = ia(29)
   l_sb = ia(30)
   l_sx = ia(31)
   l_ligw = ia(32)
   l_lrgw = ia(33)
   l_igwk = ia(34)
   l_rgwk = ia(35)
   l_nwksz = ia(36)
   l_iwksz = ia(37)
   l_rwork = ia(38)
   l_iwork = ia(39)

! Evaluate matrix-vector multiplication Ax using the fast multipole BEM

!    write (10,"(A)") "=========================== matvec: start"
!    write (10,"(/A)") "matvec, u"
!    write(10,103) (u(i),i=1,n)
! 103 format(20(ES13.6,2x))


   call upward(u, am(l_n), am(l_y), am(l_node),&
      am(l_dnorm), am(l_bc), am(l_a), am(l_xmax),&
      am(l_xmin), am(l_ymax), am(l_ymin), am(l_ielem),&
      am(l_itree), am(l_level), am(l_loct), am(l_numt),&
      am(l_ifath), am(l_nexp), am(l_ncellmx), am(l_lowlev),&
      am(l_maxl))

   call dwnwrd(u,ax, am(l_n), am(l_x), am(l_y),&
      am(l_node), am(l_dnorm), am(l_bc), am(l_a),&
      am(l_b), am(l_xmax), am(l_xmin), am(l_ymax),&
      am(l_ymin), am(l_ielem), am(l_itree), am(l_level),&
      am(l_loct), am(l_numt), am(l_ifath), am(l_nexp),&
      am(l_ntylr), am(l_ncellmx), am(l_lowlev), am(l_maxl),&
      am(l_rwork), am(l_iwork))

!    write (10,"(/A)") "matvec, ax"
!    write(10,102) (ax(i),i=1,n)
! 102 format(20(ES13.6,2x))

   return
end

!----------------------------------------------------------------------------
subroutine msolve(n,r,z,nelt,ia,ja,am,isym,rwork,iwork)
   implicit real*8(a-h,o-z)
   complex*16 am(*)
   dimension r(*),z(*),ia(*),ja(1),rwork(*),iwork(*)

! Load the pointers
   l_loct = ia(16)
   l_numt = ia(17)

! Compute the preconditioning matrix
   call msolveinv(r,z,rwork,iwork, am(l_loct), am(l_numt))
   return
end

!----------------------------------------------------------------------------
subroutine msolveinv(r,z,rwork,iwork,loct,numt)
   implicit real*8(a-h,o-z)
   dimension r(*),z(*),iwork(*),rwork(*),loct(*),numt(*)

   nleaf = iwork(1)
   do l = 1,nleaf
      l3 = l*3-1
      inod = iwork(l3)

      indr = iwork(l3+1)
      indx = loct(inod)
      indi = indx+3*nleaf+1
      nr = numt(inod)

      call dcopy(nr,r(indx),1,z(indx),1)
      call dluax(rwork(indr),nr,nr,z(indx),iwork(l3+2),&
         iwork(indi),icon)

      if(icon.ne.0) then
         write(*,*) " dluax error, icon =", icon
         stop
      endif
      iwork(l3+2) = 2
   enddo

   return
end

!----------------------------------------------------------------------------
! This subroutine solves linear system of equations Ax=b by LU decomposition.
!
! a .... given regular coefficient matrix.
! k .... given adjustable dimension for array a.
! n .... given order of matrix a.
! b .... given constant vector.
! isw .... given control information: c if 1, solve equations entirely.
! if 2, solve equations with last LU-decomposed entries.
! ip .... auxiliary 1 dimensioned array, size is n.
! transposition vector which represents
! row-exchanging by partial pivoting.
! icon.... resultant condition code.
!
! Slave subroutines used (available at www.netlib.org):
! dgetrf, dgetrs
!
!----------------------------------------------------------------------------

subroutine dluax(a,k,n,b,isw,ip,icon)
   implicit real*8 (a-h,o-z)
   dimension a(k,n),b(n),ip(n)
   data ione/1/

   icon = 30000
   if(isw.eq.1) go to 1000
   if(isw.eq.2) go to 1100
   go to 8000
1000 call dgetrf(n,n,a,k,ip,icon)
   isw = 2
1100 call dgetrs('n',n,ione,a,k,ip,b,n,icon)
8000 continue

   return
end

!----------------------------------------------------------------------------
subroutine upward(u,n,y,node,dnorm,bc,a,xmax,xmin,ymax,ymin,ielem,&
   itree,level,loct,numt,ifath,nexp,ncellmx,lowlev,maxl)

   implicit real*8(a-h,o-z)
   complex*16 a,b, z0,zi
   dimension a(0:nexp,ncellmx),&
      y(2,*),node(2,*),dnorm(2,*),bc(2,*),u(*),&
      ielem(*),itree(*),level(0:*),loct(*),numt(*),ifath(*)

   do i=1,level(lowlev+1)-1
      do k=0,nexp
         a(k,i) = (0.d0,0.d0) ! Clear multipole moments
      enddo
   enddo

   do 10 lev=lowlev,2,-1 ! Loop from leaf to level 2 cells (Upward)
      ndivx = 2**lev
      dx = (xmax-xmin)/ndivx ! Determine cell size
      dy = (ymax-ymin)/ndivx

      do 20 icell=level(lev),level(lev+1)-1 ! Loop for level l cells
         itr = itree(icell)
         itrx = mod(itr,ndivx)
         itry = itr/ndivx ! Position of the cell
         cx = xmin+(itrx + 0.5d0)*dx
         cy = ymin+(itry + 0.5d0)*dy ! Center of the cell

! Multipole expansion
         if(numt(icell).le.maxl .or. lev.eq.lowlev) then ! Compute moment
            call moment(a(0,icell),y,node,ielem(loct(icell)),&
               numt(icell),nexp,cx,cy,u(loct(icell)),&
               bc,dnorm)
         endif

! M2M translation
         if(lev.ne.2) then ! Do M2M translation to form moments
            cxp = xmin+(int(itrx/2)*2 + 1)*dx
            cyp = ymin+(int(itry/2)*2 + 1)*dy ! Center of parent cell
            z0 = cmplx(cx-cxp, cy-cyp) ! (z_c - z_c')
            io = ifath(icell) ! Cell no. of parent cell
            zi = (1.d0,0.d0)
            do k=0,nexp
               do m=k,nexp
                  a(m,io) = a(m,io) + zi*a(m-k,icell) ! Use M2M
               enddo
               zi = zi*z0/(k+1)
            enddo
         endif

20    continue
10 continue
   return
end

!----------------------------------------------------------------------------
subroutine moment(a,y,node,ielem,num,nexp,cx,cy,u,bc,dnorm)
   implicit real*8(a-h,o-y)
   implicit complex*16(z)
   complex*16 a(0:*)
   dimension y(2,*),node(2,*),ielem(*),u(*), bc(2,*),dnorm(2,*)

   do i=1,num ! Over elements in the leaf
      nelm = ielem(i) ! Element number
      n1 = node(1,nelm) ! Two ends of the element
      n2 = node(2,nelm)
      z1 = cmplx(y(1,n1)-cx, y(2,n1)-cy)
      z2 = cmplx(y(1,n2)-cx, y(2,n2)-cy)
      zwbar = conjg(z2 - z1)
      zwbar = zwbar/abs(zwbar) ! omega bar
      zp1 = z1*zwbar
      zp2 = z2*zwbar
      znorm = cmplx(dnorm(1,nelm),dnorm(2,nelm)) ! complex normal n
      if(bc(1,nelm) .eq. 1.d0) then ! Assign values to phi and q
         phi = 0.D0
         q = u(i)
      else if(bc(1,nelm) .eq. 2.d0) then
         phi = u(i)
         q = 0.D0
      endif

! Compute moments:
      a(0) = a(0) - (zp2-zp1)*q ! G kernel
      do k=1,nexp
         a(k) = a(k) + (zp2-zp1)*znorm*phi ! F kernel
         zp1 = zp1*z1/(k+1)
         zp2 = zp2*z2/(k+1)
         a(k) = a(k) - (zp2-zp1)*q ! G kernel
      enddo
   enddo
   return
end

!----------------------------------------------------------------------------
subroutine dwnwrd(u,ax,n,x,y,node,dnorm,bc,a,b,xmax,xmin,&
   ymax,ymin,ielem,itree,level,loct,numt,ifath,&
   nexp,ntylr,ncellmx,lowlev,maxl,rwork,iwork)

   implicit real*8(a-h,o-z)
   complex*16 a,b, z0,zi,zo,zp

   dimension a(0:nexp,ncellmx),b(0:ntylr,ncellmx),x(2,*),y(2,*),&
      node(2,*),dnorm(2,*),bc(2,*),u(*),ax(*),ielem(*),&
      itree(*),level(0:*),loct(*),numt(*),ifath(*),&
      rwork(*),iwork(*)

   data pi/3.141592653589793D0/
   pi2 = pi*2.d0
   do i=1,level(lowlev+1)-1
      do k=0,ntylr
         b(k,i) = (0.d0,0.d0)
      enddo
   enddo
   do i=1,n
      ax(i) = 0.d0
   enddo

   leaf = 0
   indr = 1
   indi = 1
   do 110 lev=2,lowlev ! Downward from level 2 cells to leaf cells
      ndivx = 2**lev
      dx = (xmax-xmin)/ndivx
      dy = (ymax-ymin)/ndivx
      do 120 icell=level(lev),level(lev+1)-1 ! Loop for level l cells
         itr = itree(icell)
         itrx = mod(itr,ndivx)
         itry = itr/ndivx ! Position of the cell
         cx = xmin+(itrx + 0.5d0)*dx
         cy = ymin+(itry + 0.5d0)*dy ! Center of the cell
         itrxp = itrx/2
         itryp = itry/2

! From the parent cell (use L2L)
         if(lev.ne.2) then
            cxp = xmin+(itrxp*2+1)*dx
            cyp = ymin+(itryp*2+1)*dy ! Center of the parent cell
            z0 = cmplx(cx-cxp, cy-cyp)
            io = ifath(icell) ! Cell no. of the parent cell
            zi = (1.d0,0.d0)
            do k=0,ntylr
               do m=0,ntylr-k

                  b(m,icell) = b(m,icell) + zi*b(k+m,io) ! L2L translation
               enddo
               zi = zi*z0/(k+1)
            enddo
         endif

         do 130 jcell=level(lev),level(lev+1)-1
            jtr = itree(jcell)
            jtrx = mod(jtr,ndivx)
            jtry = jtr/ndivx
            jtrxp = jtrx/2
            jtryp = jtry/2

! The parents must be neighbours
            if(iabs(itrxp-jtrxp).gt.1 .or. iabs(itryp-jtryp).gt.1)&
               goto 130

! For non-neighbours (cells in interaction list) (use M2L)
            if(iabs(itrx-jtrx).gt.1 .or. iabs(itry-jtry).gt.1) then
               ccx = xmin + (jtrx + 0.5d0)*dx
               ccy = ymin + (jtry + 0.5d0)*dy ! Center of the j cell
               z0 = cmplx(cx-ccx, cy-ccy)
               b(0,icell) = b(0,icell) - log(z0)*a(0,jcell)
               zo = 1.
               do m=1,nexp+ntylr
                  zo = zo/z0
                  kmin = max(0,m-nexp)
                  kmax = min(m,ntylr)
                  sgn = (-1.0)**kmin
                  do k=kmin,kmax
                     b(k,icell) = b(k,icell) + sgn*zo*a(m-k,jcell) ! M2L
                     sgn = -sgn
                  enddo
                  zo = zo*m
               enddo

! Contribution from neighbouring leaves (use direct)
            elseif(numt(jcell).le.maxl .or. &
               numt(icell).le.maxl .or. lev.eq.lowlev) then
               if(icell.eq.jcell) then
                  leaf = leaf+1
                  leaf3 = leaf*3-1
                  if(iwork(leaf3).ne.icell) then
                     write(3,*) leaf,iwork(1),iwork(leaf3),icell,'check'
                     write(*,*) " icell error"
                     stop
                  endif
                  indr = iwork(leaf3+1)
                  indi = iwork(leaf3+2)
               endif

               call direct(ielem(loct(icell)),ielem(loct(jcell)),&
                  node,x,y,numt(icell),numt(jcell),dnorm,&
                  ax(loct(icell)),u(loct(jcell)),icell,jcell,&
                  rwork(indr),indi,bc) ! Direct integration
            endif
130      continue

! Compute Ax if reach a leaf (Evaluate local expansion at each
! collocation point)
         if(numt(icell).le.maxl .or. lev.eq.lowlev) then
            fact = 1.d0
            do itylr=1,ntylr
               fact = fact/itylr
               b(itylr,icell) = b(itylr,icell)*fact
            enddo

            do in=1,numt(icell)
               inax = loct(icell) + in-1 ! Element number in the tree
               indx = ielem(inax) ! Original element number
               zp = b(ntylr,icell)
               z0 = cmplx(x(1,indx)-cx, x(2,indx)-cy)
               do itylr=ntylr-1,0,-1
                  zp = zp*z0 + b(itylr,icell) ! Local expansion
               enddo
               zp = zp/pi2
               ax(inax) = ax(inax) + dreal(zp) ! Array Ax
            enddo
         endif

120   continue
110 continue

   return
end

!----------------------------------------------------------------------------
subroutine direct(inod,jnod,node,x,y,ni,nj,dnorm,ax,u,icell,jcell, amat,isw,bc)
   implicit real*8(a-h,o-z)
   dimension inod(*),jnod(*),node(2,*),x(2,*),y(2,*),&
      dnorm(2,*),ax(*),u(*),amat(ni,*), bc(2,*)
   data pi/3.141592653589793D0/
   pi2 = pi*2.d0
   do j = 1, nj

      jind = jnod(j)
      al = sqrt((y(1,node(1,jind))- y(1,node(2,jind)))**2 +&
         (y(2,node(1,jind)) - y(2,node(2,jind)))**2) ! Element length

      do i = 1, ni
         iind = inod(i)
         x11 = y(1,node(1,jind)) - x(1,iind)
         x21 = y(2,node(1,jind)) - x(2,iind)
         x12 = y(1,node(2,jind)) - x(1,iind)
         x22 = y(2,node(2,jind)) - x(2,iind)
         r1 = sqrt(x11**2 + x21**2)
         r2 = sqrt(x12**2 + x22**2)
         d = x11*dnorm(1,jind) + x21*dnorm(2,jind)
         t1 = -x11*dnorm(2,jind) + x21*dnorm(1,jind)
         t2 = -x12*dnorm(2,jind) + x22*dnorm(1,jind)
         ds = abs(d)
         dtheta = datan2(ds*al,ds**2+t1*t2)
         aa = (-dtheta*ds + al + t1*log(r1) - t2*log(r2))/pi2
         if(d .lt. 0.d0) dtheta = -dtheta

         bb = -dtheta/pi2
         if(iind .eq. jind) bb = 0.5d0
         if(bc(1,jind) .eq. 1.) then
            ax(i) = ax(i) - aa*u(j)
            if(icell.eq.jcell .and. isw.eq.1) then ! Store coefficients in
               amat(i,j) = - aa ! first iteration
            endif
         else if(bc(1,jind) .eq. 2.) then
            ax(i) = ax(i) + bb*u(j)
            if(icell.eq.jcell .and. isw.eq.1) then ! Store coefficients in
               amat(i,j) = bb ! first iteration
            endif
         endif
      enddo
   enddo
   return
end

!----------------------------------------------------------------------------
subroutine domain_field(nfield,xfield,f,n,x,y,bc,node,dnorm,u)
   implicit real*8(a-h,o-z)

   dimension xfield(2,*), f(*), x(2,*),y(2,*),bc(2,*),node(2,*),&
      dnorm(2,*),u(*)
   data pi/3.141592653589793D0/
   pi2 = pi*2.d0
   do i=1,nfield
      f(i) = 0.d0
   enddo
   do j=1,n ! Loop over all elements
      if(bc(1,j).eq.1) then
         f0 = bc(2,j)
         df0 = u(j)

      else if(bc(1,j).eq.2) then

         f0 = u(j)
         df0 = bc(2,j)
      endif
      al= sqrt((y(1,node(2,j)) - y(1,node(1,j)))**2 +&
         (y(2,node(2,j)) - y(2,node(1,j)))**2) ! Element length

      do i=1,nfield ! Loop over all field points inside the domain
         x11 = y(1,node(1,j)) - xfield(1,i)
         x21 = y(2,node(1,j)) - xfield(2,i)
         x12 = y(1,node(2,j)) - xfield(1,i)
         x22 = y(2,node(2,j)) - xfield(2,i)
         r1 = sqrt(x11**2 + x21**2)
         r2 = sqrt(x12**2 + x22**2)
         d = x11*dnorm(1,j) + x21*dnorm(2,j)
         t1 = -x11*dnorm(2,j) + x21*dnorm(1,j)
         t2 = -x12*dnorm(2,j) + x22*dnorm(1,j)
         ds = abs(d)
         theta1 = datan2(t1,ds)
         theta2 = datan2(t2,ds)
         dtheta = theta2 - theta1

         aa = (-dtheta*ds + al + t1*log(r1) - t2*log(r2))/pi2
         if(d .lt. 0.d0) dtheta = -dtheta
         bb = -dtheta/pi2
         f(i) = f(i) + aa*df0 - bb*f0
      enddo
   enddo

! Output results
   do i=1,nfield
      write(9,20) xfield(1,i), f(i)
   enddo
20 format(1x, 4E18.8)
   return
end

!----------------------------------------------------------------------------
